From: Drew Parsons Date: Wed, 14 Nov 2018 15:10:51 +0000 (+0800) Subject: Import pygalmesh_0.2.6.orig.tar.gz X-Git-Tag: archive/raspbian/0.10.6-5+rpi1~1^2^2^2~1 X-Git-Url: https://dgit.raspbian.org/%22http://www.example.com/cgi/success/%22http:/www.example.com/cgi/success?a=commitdiff_plain;h=0a97edc2cc0a5b9a27c9302c59b5273c7179b4ab;p=pygalmesh.git Import pygalmesh_0.2.6.orig.tar.gz [dgit import orig pygalmesh_0.2.6.orig.tar.gz] --- 0a97edc2cc0a5b9a27c9302c59b5273c7179b4ab diff --git a/.circleci/config.yml b/.circleci/config.yml new file mode 100644 index 0000000..de4f763 --- /dev/null +++ b/.circleci/config.yml @@ -0,0 +1,21 @@ +version: 2 + +jobs: + + build: + working_directory: ~/pygalmesh + docker: + - image: ubuntu:18.04 + steps: + - run: apt-get update + - run: apt-get install -y libcgal-dev libeigen3-dev python3-pip clang++-6.0 + - run: pip3 install -U pytest pytest-cov black flake8 meshio pybind11 + - checkout + # install + # Use clang++ for its smaller memory footprint. + - run: CC="clang++-6.0" pip3 install . + # lint + - run: LC_ALL=C.UTF-8 black --check setup.py pygalmesh/ test/*.py + - run: flake8 setup.py pygalmesh/ test/*.py + # The actual test + - run: cd test/ && pytest diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..c321e71 --- /dev/null +++ b/.flake8 @@ -0,0 +1,5 @@ +[flake8] +ignore = E203, E266, E501, W503 +max-line-length = 80 +max-complexity = 18 +select = B,C,E,F,W,T4,B9 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..62bde14 --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +*.off +*.vtu +*.mesh +.cache/ +build/ +dist/ +MANIFEST +README.rst +do-configure.sh +.pytest_cache/ +*.egg-info/ diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 0000000..31f545d --- /dev/null +++ b/.pylintrc @@ -0,0 +1,8 @@ +[MESSAGES CONTROL] + +disable= + fixme, + invalid-name, + missing-docstring, + no-member, + too-few-public-methods diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..69ee8c6 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,34 @@ +dist: trusty + +language: python + +python: + - '2.7' + - '3.6' + +addons: + apt: + sources: + - sourceline: 'ppa:nschloe/cgal-backports' + - sourceline: 'ppa:nschloe/eigen-backports' + - sourceline: 'ppa:nschloe/swig-backports' + packages: + - libcgal-dev + - libcgal-qt5-11 # bug in cgal, should be installed automatically + - libeigen3-dev + - pandoc + - python-numpy + - swig + +# test dependencies +before_install: + - pip install meshio + +install: + - pip install pytest + - pip install . + +script: + - cd test + - curl https://raw.githubusercontent.com/libigl/libigl-unit-tests/master/data/elephant.off > elephant.off + - pytest diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..fb0710f --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required(VERSION 3.0) + +project(pygalmesh CXX) + +add_subdirectory(src) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..281c7b5 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2016-2018 Nico Schlömer + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..c2fb0da --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,6 @@ +include src/domain.hpp +include src/generate_from_off.hpp +include src/generate.hpp +include src/generate_surface_mesh.hpp +include src/polygon2d.hpp +include src/primitives.hpp diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..57e7287 --- /dev/null +++ b/Makefile @@ -0,0 +1,32 @@ +VERSION=$(shell python3 -c "import pygalmesh; print(pygalmesh.__version__)") + +default: + @echo "\"make publish\"?" + +tag: + @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi + @echo "Tagging v$(VERSION)..." + git tag v$(VERSION) + git push --tags + +upload: setup.py + # Make sure we're on the master branch + @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi + rm -rf dist/* + python3 setup.py sdist + twine upload dist/*.tar.gz + # HTTPError: 400 Client Error: Binary wheel 'pygalmesh-0.2.0-cp27-cp27mu-linux_x86_64.whl' has an unsupported platform tag 'linux_x86_64'. for url: https://upload.pypi.org/legacy/ + # python3 setup.py bdist_wheel --universal + # twine upload dist/*.whl + +publish: tag upload + +clean: + @find . | grep -E "(__pycache__|\.pyc|\.pyo$\)" | xargs rm -rf + +black: + black setup.py pygalmesh/ test/*.py + +lint: + black --check setup.py pygalmesh/ test/*.py + flake8 setup.py pygalmesh/ test/*.py diff --git a/README.md b/README.md new file mode 100644 index 0000000..fc45727 --- /dev/null +++ b/README.md @@ -0,0 +1,321 @@ +# pygalmesh + +A Python frontend to [CGAL](https://www.cgal.org/)'s [3D mesh generation +capabilities](https://doc.cgal.org/latest/Mesh_3/index.html). + +[![CircleCI](https://img.shields.io/circleci/project/github/nschloe/pygalmesh/master.svg)](https://circleci.com/gh/nschloe/pygalmesh/tree/master) +[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black) +[![PyPi Version](https://img.shields.io/pypi/v/pygalmesh.svg)](https://pypi.org/project/pygalmesh) +[![GitHub stars](https://img.shields.io/github/stars/nschloe/pygalmesh.svg?label=Stars&logo=github)](https://github.com/nschloe/pygalmesh) + +pygalmesh makes it easy to create high-quality 3D volume and surface meshes. + +### Background + +CGAL offers two different approaches for mesh generation: + +1. Meshes defined implicitly by level sets of functions. +2. Meshes defined by a set of bounding planes. + +pygalmesh provides a front-end to the first approach, which has the following +advantages and disadvantages: + +* All boundary points are guaranteed to be in the level set within any specified + residual. This results in smooth curved surfaces. +* Sharp intersections of subdomains (e.g., in unions or differences of sets) + need to be specified manually (via feature edges, see below), which can be + tedious. + +On the other hand, the bounding-plane approach (realized by +[mshr](https://bitbucket.org/fenics-project/mshr)), has the following +properties: + +* Smooth, curved domains are approximated by a set of bounding planes, + resulting in more of less visible edges. +* Intersections of domains can be computed automatically, so domain unions etc. + have sharp edges where they belong. + +Other Python mesh generators are [pygmsh](https://github.com/nschloe/pygmsh) (a +frontend to [gmsh](http://gmsh.info/)) and +[MeshPy](https://github.com/inducer/meshpy). +[meshzoo](https://github.com/nschloe/meshzoo) provides some basic canonical +meshes. + +### Examples + +#### A simple ball + + +```python +import pygalmesh + +s = pygalmesh.Ball([0, 0, 0], 1.0) +pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2) +``` +CGAL's mesh generator returns Medit-files, which can be processed by, e.g., +[meshio](https://github.com/nschloe/meshio). +```python +import meshio +vertices, cells, _, _, _ = meshio.read("out.mesh") +``` +The mesh generation comes with many more options, described +[here](https://doc.cgal.org/latest/Mesh_3/). Try, for example, +```python +pygalmesh.generate_mesh( + s, + "out.mesh", + cell_size=0.2, + edge_size=0.1, + odt=True, + lloyd=True, + verbose=False +) +``` + +#### Other primitive shapes + + +pygalmesh provides out-of-the-box support for balls, cuboids, ellipsoids, tori, +cones, cylinders, and tetrahedra. Try for example +```python +import pygalmesh + +s0 = pygalmesh.Tetrahedron( + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0] +) +pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1) +``` + +#### Domain combinations + + +Supported are unions, intersections, and differences of all domains. As +mentioned above, however, the sharp intersections between two domains are not +automatically handled. Try for example +```python +import pygalmesh + +radius = 1.0 +displacement = 0.5 +s0 = pygalmesh.Ball([displacement, 0, 0], radius) +s1 = pygalmesh.Ball([-displacement, 0, 0], radius) +u = pygalmesh.Difference(s0, s1) +``` +To sharpen the intersection circle, add it as a feature edge polygon line, +e.g., +```python +a = numpy.sqrt(radius**2 - displacement**2) +edge_size = 0.15 +n = int(2*numpy.pi*a / edge_size) +circ = [ + [ + 0.0, + a * numpy.cos(i * 2*numpy.pi / n), + a * numpy.sin(i * 2*numpy.pi / n) + ] for i in range(n) + ] +circ.append(circ[0]) + +pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + facet_angle=25, + facet_size=0.15, + cell_radius_edge_ratio=2.0 +) +``` +Note that the length of the polygon legs are kept in sync with the `edge_size` +of the mesh generation. This makes sure that it fits in nicely with the rest of +the mesh. + +#### Domain deformations + + +You can of course translate, rotate, scale, and stretch any domain. Try, for +example, +```python +import pygalmesh + +s = pygalmesh.Stretch( + pygalmesh.Ball([0, 0, 0], 1.0), + [1.0, 2.0, 0.0] +) + +pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.1) +``` + +#### Extrusion of 2D polygons + + +pygalmesh lets you extrude any polygon into a 3D body. It even supports +rotation alongside! +```python +import pygalmesh + +p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) +edge_size = 0.1 +domain = pygalmesh.Extrude( + p, + [0.0, 0.0, 1.0], + 0.5 * 3.14159265359, + edge_size +) +pygalmesh.generate_mesh( + domain, + "out.mesh", + cell_size=0.1, + edge_size=edge_size, + verbose=False +) +``` +Feature edges are automatically preserved here, which is why an edge length +needs to be given to `pygalmesh.Extrude`. + +#### Rotation bodies + + +Polygons in the x-z-plane can also be rotated around the z-axis to yield a +rotation body. +```python +import pygalmesh + +p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]]) +edge_size = 0.1 +domain = pygalmesh.RingExtrude(p, edge_size) +pygalmesh.generate_mesh( + domain, + "out.mesh", + cell_size=0.1, + edge_size=edge_size, + verbose=False +) +``` + +#### Your own custom level set function + + +If all of the variety is not enough for you, you can define your own custom +level set function. You simply need to subclass `pygalmesh.DomainBase` and +specify a function, e.g., +```python +import pygalmesh +class Heart(pygalmesh.DomainBase): + def __init__(self): + super(Heart, self).__init__() + return + + def eval(self, x): + return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \ + - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3 + + def get_bounding_sphere_squared_radius(self): + return 10.0 + +d = Heart() +pygalmesh.generate_mesh(d, "out.mesh", cell_size=0.1) +``` +Note that you need to specify the square of a bounding sphere radius, used as +an input to CGAL's mesh generator. + +#### Surface meshes + +If you're only after the surface of a body, pygalmesh has +`generate_surface_mesh` for you. It offers fewer options (obviously, +`cell_size` is gone), but otherwise works the same way: +```python +import pygalmesh + +s = pygalmesh.Ball([0, 0, 0], 1.0) +pygalmesh.generate_surface_mesh( + s, + "out.off", + angle_bound=30, + radius_bound=0.1, + distance_bound=0.1 +) +``` +The output format is +[OFF](http://segeval.cs.princeton.edu/public/off_format.html) which again is +handled by [meshio](https://github.com/nschloe/meshio). + +Refer to [CGAL's +documention](https://doc.cgal.org/latest/Surface_mesher/index.html) for the +options. + +#### Meshes from OFF files + + +If you have an OFF file at hand (like +[elephant.off](https://raw.githubusercontent.com/CGAL/cgal-swig-bindings/master/examples/data/elephant.off) +or [these](https://github.com/CGAL/cgal/tree/master/Surface_mesher/demo/Surface_mesher/inputs)), +pygalmesh generates the mesh via +```python +import pygalmesh + +pygalmesh.generate_from_off( + "elephant.off", + "out.mesh", + facet_angle=25.0, + facet_size=0.15, + facet_distance=0.008, + cell_radius_edge_ratio=3.0, + verbose=False +) +``` + +### Installation + +For installation, pygalmesh needs [CGAL](https://www.cgal.org/) and +[Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) installed on your +system. They are typically available on your Linux distribution, e.g., on +Ubuntu +``` +sudo apt install libcgal-dev libeigen3-dev +``` +After that, pygalmesh can be [installed from the Python Package +Index](https://pypi.org/project/pygalmesh/), so with +``` +pip install -U pygalmesh +``` +you can install/upgrade. + +[meshio](https://github.com/nschloe/meshio) (`sudo -H pip install meshio`) +can be helpful in processing the meshes. + +#### Manual installation + +For manual installation (if you're a developer or just really keen on getting +the bleeding edge version of pygalmesh), there are two possibilities: + + * Get the sources, type `sudo python setup.py install`. This does the trick + most the time. + * As a fallback, there's a CMake-based installation. Simply go `cmake + /path/to/sources/` and `make`. + +### Testing + +To run the pygalmesh unit tests, check out this repository and type +``` +pytest +``` + +### Distribution + +To create a new release + +1. bump the `__version__` number (in `setup.py` _and_ `src/pygalmesh.i`) + +2. publish to PyPi and GitHub: + ``` + make publish + ``` + +### License + +pygalmesh is published under the [MIT license](https://en.wikipedia.org/wiki/MIT_License). diff --git a/pygalmesh/__about__.py b/pygalmesh/__about__.py new file mode 100644 index 0000000..5bf2aa0 --- /dev/null +++ b/pygalmesh/__about__.py @@ -0,0 +1,11 @@ +# -*- coding: utf-8 -*- +# + +__author__ = u"Nico Schlömer" +__author_email__ = "nico.schloemer@gmail.com" +__copyright__ = u"Copyright (c) 2016-2018, {} <{}>".format(__author__, __author_email__) +__license__ = "License :: OSI Approved :: MIT License" +__version__ = "0.2.6" +__maintainer__ = u"Nico Schlömer" +__status__ = "Development Status :: 4 - Beta" +__url__ = "https://github.com/nschloe/pygalmesh" diff --git a/pygalmesh/__init__.py b/pygalmesh/__init__.py new file mode 100644 index 0000000..dc00c34 --- /dev/null +++ b/pygalmesh/__init__.py @@ -0,0 +1,70 @@ +# -*- coding: utf-8 -*- +# +# https://github.com/pybind/pybind11/issues/1004 +from _pygalmesh import ( + DomainBase, + Translate, + Rotate, + Scale, + Stretch, + Intersection, + Union, + Difference, + Extrude, + Ball, + Cuboid, + Ellipsoid, + Tetrahedron, + Cone, + Cylinder, + Torus, + HalfSpace, + Polygon2D, + RingExtrude, + generate_mesh, + generate_surface_mesh, + generate_from_off, +) + +from .__about__ import ( + __author__, + __author_email__, + __copyright__, + __license__, + __version__, + __maintainer__, + __status__, +) + +__all__ = [ + "__author__", + "__author_email__", + "__copyright__", + "__license__", + "__version__", + "__maintainer__", + "__status__", + # + "DomainBase", + "Translate", + "Rotate", + "Scale", + "Stretch", + "Intersection", + "Union", + "Difference", + "Extrude", + "Ball", + "Cuboid", + "Ellipsoid", + "Tetrahedron", + "Cone", + "Cylinder", + "Torus", + "HalfSpace", + "Polygon2D", + "RingExtrude", + "generate_mesh", + "generate_surface_mesh", + "generate_from_off", +] diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..16663a8 --- /dev/null +++ b/setup.py @@ -0,0 +1,84 @@ +# -*- coding: utf-8 -*- +# +import os + +import codecs +from setuptools import setup, Extension, find_packages + + +# https://packaging.python.org/single_source_version/ +base_dir = os.path.abspath(os.path.dirname(__file__)) +about = {} +with open(os.path.join(base_dir, "pygalmesh", "__about__.py"), "rb") as handle: + # pylint: disable=exec-used + exec(handle.read(), about) + + +class get_pybind_include(object): + """Helper class to determine the pybind11 include path + The purpose of this class is to postpone importing pybind11 + until it is actually installed, so that the ``get_include()`` + method can be invoked. + """ + + def __init__(self, user=False): + self.user = user + + def __str__(self): + import pybind11 + + return pybind11.get_include(self.user) + + +def read(fname): + return codecs.open(os.path.join(base_dir, fname), encoding="utf-8").read() + + +ext_modules = [ + Extension( + "_pygalmesh", + [ + "src/generate.cpp", + "src/generate_from_off.cpp", + "src/generate_surface_mesh.cpp", + "src/pybind11.cpp", + ], + language="c++", + include_dirs=[ + "/usr/include/eigen3/", + # Path to pybind11 headers + get_pybind_include(), + get_pybind_include(user=True), + ], + libraries=["CGAL", "gmp", "mpfr"], + # extra_compile_args=['-std=c++11'] + ) +] + +setup( + name="pygalmesh", + packages=find_packages(), + # cmdclass={'build_ext': BuildExt}, + ext_modules=ext_modules, + version=about["__version__"], + url=about["__url__"], + author=about["__author__"], + author_email=about["__author_email__"], + install_requires=["numpy", "pybind11 >= 2.2", "pipdate"], + description="Python frontend to CGAL's 3D mesh generation capabilities", + long_description=read("README.md"), + long_description_content_type="text/markdown", + license=about["__license__"], + classifiers=[ + about["__status__"], + about["__license__"], + "Operating System :: OS Independent", + "Programming Language :: Python", + "Programming Language :: Python :: 2", + "Programming Language :: Python :: 3", + "Topic :: Scientific/Engineering", + "Topic :: Scientific/Engineering :: Mathematics", + "Topic :: Scientific/Engineering :: Physics", + "Topic :: Scientific/Engineering :: Visualization", + ], +) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..a6df64b --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,31 @@ +FIND_PACKAGE(pybind11 REQUIRED) +# include_directories(${PYBIND11_INCLUDE_DIR}) + +FIND_PACKAGE(Eigen3 REQUIRED) +include_directories(${EIGEN3_INCLUDE_DIR}) + +FILE(GLOB pygalmesh_SRCS "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp") +# FILE(GLOB pygalmesh_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp") + +pybind11_add_module(pygalmesh ${pygalmesh_SRCS}) + +# Call CGAL as late as possible. It messes around with some +# compiler variables. See bugs +# and +# . +FIND_PACKAGE(CGAL REQUIRED) +include(${CGAL_USE_FILE}) + +# ADD_LIBRARY(pygalmesh ${pygalmesh_SRCS}) +target_link_libraries(pygalmesh PRIVATE ${CGAL_LIBRARIES}) + +# execute_process( +# COMMAND python -c "from distutils.sysconfig import get_python_lib; print get_python_lib()" +# OUTPUT_VARIABLE PYTHON_SITE_PACKAGES +# OUTPUT_STRIP_TRAILING_WHITESPACE +# ) +# install(TARGETS _pygalmesh DESTINATION ${PYTHON_SITE_PACKAGES}) +# install( +# FILES ${CMAKE_BINARY_DIR}/src/pygalmesh.py +# DESTINATION ${PYTHON_SITE_PACKAGES} +# ) diff --git a/src/domain.hpp b/src/domain.hpp new file mode 100644 index 0000000..ff8cd9c --- /dev/null +++ b/src/domain.hpp @@ -0,0 +1,490 @@ +#ifndef DOMAIN_HPP +#define DOMAIN_HPP + +#include +#include +#include +#include +#include + +namespace pygalmesh { + +class DomainBase +{ + public: + + virtual ~DomainBase() = default; + + virtual + double + eval(const std::array & x) const = 0; + + virtual + double + get_bounding_sphere_squared_radius() const = 0; + + virtual + std::vector>> + get_features() const + { + return {}; + }; +}; + +class Translate: public pygalmesh::DomainBase +{ + public: + Translate( + const std::shared_ptr & domain, + const std::array & direction + ): + domain_(domain), + direction_(Eigen::Vector3d(direction.data())), + translated_features_(translate_features(domain->get_features(), direction_)) + { + } + + virtual ~Translate() = default; + + std::vector>> + translate_features( + const std::vector>> & features, + const Eigen::Vector3d & direction + ) const + { + std::vector>> translated_features; + for (const auto & feature: features) { + std::vector> translated_feature; + for (const auto & point: feature) { + const std::array translated_point = { + point[0] + direction[0], + point[1] + direction[1], + point[2] + direction[2] + }; + translated_feature.push_back(translated_point); + } + translated_features.push_back(translated_feature); + } + return translated_features; + } + + virtual + double + eval(const std::array & x) const + { + const std::array d = { + x[0] - direction_[0], + x[1] - direction_[1], + x[2] - direction_[2] + }; + return domain_->eval(d); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double radius = sqrt(domain_->get_bounding_sphere_squared_radius()); + const double dir_norm = direction_.norm(); + return (radius + dir_norm)*(radius + dir_norm); + } + + virtual + std::vector>> + get_features() const + { + return translated_features_; + }; + + private: + const std::shared_ptr domain_; + const Eigen::Vector3d direction_; + const std::vector>> translated_features_; +}; + +class Rotate: public pygalmesh::DomainBase +{ + public: + Rotate( + const std::shared_ptr & domain, + const std::array & axis, + const double angle + ): + domain_(domain), + normalized_axis_(Eigen::Vector3d(axis.data()).normalized()), + sinAngle_(sin(angle)), + cosAngle_(cos(angle)), + rotated_features_(rotate_features(domain_->get_features())) + { + } + + virtual ~Rotate() = default; + + Eigen::Vector3d + rotate( + const Eigen::Vector3d & vec, + const Eigen::Vector3d & axis, + const double sinAngle, + const double cosAngle + ) const + { + // Rotate a vector `v` by the angle `theta` in the plane perpendicular + // to the axis given by `u`. + // Refer to + // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle + // + // cos(theta) * I * v + // + sin(theta) u\cross v + // + (1-cos(theta)) (u*u^T) * v + return cosAngle * vec + + sinAngle * axis.cross(vec) + + (1.0-cosAngle) * axis.dot(vec) * axis; + } + + virtual + double + eval(const std::array & x) const + { + // rotate with negative angle + const auto p2 = rotate( + Eigen::Vector3d(x.data()), + normalized_axis_, + -sinAngle_, + cosAngle_ + ); + return domain_->eval({p2[0], p2[1], p2[2]}); + } + + std::vector>> + rotate_features( + const std::vector>> & features + ) const + { + std::vector>> rotated_features; + for (const auto & feature: features) { + std::vector> rotated_feature; + for (const auto & point: feature) { + const auto p2 = rotate( + Eigen::Vector3d(point.data()), + normalized_axis_, + sinAngle_, + cosAngle_ + ); + rotated_feature.push_back({p2[0], p2[1], p2[2]}); + } + rotated_features.push_back(rotated_feature); + } + return rotated_features; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return domain_->get_bounding_sphere_squared_radius(); + } + + virtual + std::vector>> + get_features() const + { + return rotated_features_; + }; + + private: + const std::shared_ptr domain_; + const Eigen::Vector3d normalized_axis_; + const double sinAngle_; + const double cosAngle_; + const std::vector>> rotated_features_; +}; + +class Scale: public pygalmesh::DomainBase +{ + public: + Scale( + std::shared_ptr & domain, + const double alpha + ): + domain_(domain), + alpha_(alpha), + scaled_features_(scale_features(domain_->get_features())) + { + assert(alpha_ > 0.0); + } + + virtual ~Scale() = default; + + virtual + double + eval(const std::array & x) const + { + return domain_->eval({x[0]/alpha_, x[1]/alpha_, x[2]/alpha_}); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius(); + } + + std::vector>> + scale_features( + const std::vector>> & features + ) const + { + std::vector>> scaled_features; + for (const auto & feature: features) { + std::vector> scaled_feature; + for (const auto & point: feature) { + scaled_feature.push_back({ + alpha_ * point[0], + alpha_ * point[1], + alpha_ * point[2] + }); + } + scaled_features.push_back(scaled_feature); + } + return scaled_features; + } + + virtual + std::vector>> + get_features() const + { + return scaled_features_; + }; + + private: + std::shared_ptr domain_; + const double alpha_; + const std::vector>> scaled_features_; +}; + +class Stretch: public pygalmesh::DomainBase +{ + public: + Stretch( + std::shared_ptr & domain, + const std::array & direction + ): + domain_(domain), + normalized_direction_(Eigen::Vector3d(direction.data()).normalized()), + alpha_(Eigen::Vector3d(direction.data()).norm()), + stretched_features_(stretch_features(domain_->get_features())) + { + assert(alpha_ > 0.0); + } + + virtual ~Stretch() = default; + + virtual + double + eval(const std::array & x) const + { + const Eigen::Vector3d v(x.data()); + const double beta = normalized_direction_.dot(v); + // scale the component of normalized_direction_ by 1/alpha_ + const auto v2 = beta/alpha_ * normalized_direction_ + + (v - beta * normalized_direction_); + return domain_->eval({v2[0], v2[1], v2[2]}); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius(); + } + + std::vector>> + stretch_features( + const std::vector>> & features + ) const + { + std::vector>> stretched_features; + for (const auto & feature: features) { + std::vector> stretched_feature; + for (const auto & point: feature) { + // scale the component of normalized_direction_ by alpha_ + const Eigen::Vector3d v(point.data()); + const double beta = normalized_direction_.dot(v); + const auto v2 = beta * alpha_ * normalized_direction_ + + (v - beta * normalized_direction_); + stretched_feature.push_back({v2[0], v2[1], v2[2]}); + } + stretched_features.push_back(stretched_feature); + } + return stretched_features; + } + + virtual + std::vector>> + get_features() const + { + return stretched_features_; + }; + + private: + std::shared_ptr domain_; + const Eigen::Vector3d normalized_direction_; + const double alpha_; + const std::vector>> stretched_features_; +}; + +class Intersection: public pygalmesh::DomainBase +{ + public: + explicit Intersection( + std::vector> & domains + ): + domains_(domains) + { + } + + virtual ~Intersection() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a differentiable expression + double maxval = std::numeric_limits::lowest(); + for (const auto & domain: domains_) { + maxval = std::max(maxval, domain->eval(x)); + } + return maxval; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double min = std::numeric_limits::max(); + for (const auto & domain: domains_) { + min = std::min(min, domain->get_bounding_sphere_squared_radius()); + } + return min; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + for (const auto & domain: domains_) { + const auto f = domain->get_features(); + features.insert(std::end(features), std::begin(f), std::end(f)); + } + return features; + }; + + private: + std::vector> domains_; +}; + +class Union: public pygalmesh::DomainBase +{ + public: + explicit Union( + std::vector> & domains + ): + domains_(domains) + { + } + + virtual ~Union() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a differentiable expression + double minval = std::numeric_limits::max(); + for (const auto & domain: domains_) { + minval = std::min(minval, domain->eval(x)); + } + return minval; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & domain: domains_) { + max = std::max(max, domain->get_bounding_sphere_squared_radius()); + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + for (const auto & domain: domains_) { + const auto f = domain->get_features(); + features.insert(std::end(features), std::begin(f), std::end(f)); + } + return features; + }; + + private: + std::vector> domains_; +}; + +class Difference: public pygalmesh::DomainBase +{ + public: + Difference( + std::shared_ptr & domain0, + std::shared_ptr & domain1 + ): + domain0_(domain0), + domain1_(domain1) + { + } + + virtual ~Difference() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a continuous (perhaps even differentiable) expression + const double val0 = domain0_->eval(x); + const double val1 = domain1_->eval(x); + return (val0 < 0.0 && val1 >= 0.0) ? val0 : std::max(val0, -val1); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return domain0_->get_bounding_sphere_squared_radius(); + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + + const auto f0 = domain0_->get_features(); + features.insert(std::end(features), std::begin(f0), std::end(f0)); + + const auto f1 = domain1_->get_features(); + features.insert(std::end(features), std::begin(f1), std::end(f1)); + + return features; + }; + + private: + std::shared_ptr domain0_; + std::shared_ptr domain1_; +}; + +} // namespace pygalmesh +#endif // DOMAIN_HPP diff --git a/src/generate.cpp b/src/generate.cpp new file mode 100644 index 0000000..0e2e45f --- /dev/null +++ b/src/generate.cpp @@ -0,0 +1,144 @@ +#define CGAL_MESH_3_VERBOSE 1 + +#include "generate.hpp" + +#include + +#include +#include +#include + +#include +#include +#include + +namespace pygalmesh { + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +// Wrapper for DomainBase for translating to K::Point_3/FT. +class CgalDomainWrapper +{ + public: + explicit CgalDomainWrapper(const std::shared_ptr & domain): + domain_(domain) + { + } + + virtual ~CgalDomainWrapper() = default; + + virtual + K::FT + operator()(K::Point_3 p) const + { + return domain_->eval({p.x(), p.y(), p.z()}); + } + + private: + const std::shared_ptr domain_; +}; + +typedef CGAL::Mesh_domain_with_polyline_features_3< + CGAL::Implicit_mesh_domain_3 + > + Mesh_domain; + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Mesh Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; +typedef Mesh_criteria::Facet_criteria Facet_criteria; +typedef Mesh_criteria::Cell_criteria Cell_criteria; + +// translate vector> to list> +std::list> +translate_feature_edges( + const std::vector>> & feature_edges + ) +{ + std::list> polylines; + for (const auto & feature_edge: feature_edges) { + std::vector polyline; + for (const auto & point: feature_edge) { + polyline.push_back(K::Point_3(point[0], point[1], point[2])); + } + polylines.push_back(polyline); + } + return polylines; +} + +void +generate_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const std::vector>> & feature_edges, + const double bounding_sphere_radius, + const double boundary_precision, + const bool lloyd, + const bool odt, + const bool perturb, + const bool exude, + const double edge_size, + const double facet_angle, + const double facet_size, + const double facet_distance, + const double cell_radius_edge_ratio, + const double cell_size, + const bool verbose + ) +{ + const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ? + bounding_sphere_radius*bounding_sphere_radius : + // some wiggle room + 1.01 * domain->get_bounding_sphere_squared_radius(); + + const auto d = CgalDomainWrapper(domain); + Mesh_domain cgal_domain( + d, + K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2), + boundary_precision + ); + + const auto native_features = translate_feature_edges(domain->get_features()); + cgal_domain.add_features(native_features.begin(), native_features.end()); + + const auto polylines = translate_feature_edges(feature_edges); + cgal_domain.add_features(polylines.begin(), polylines.end()); + + Mesh_criteria criteria( + CGAL::parameters::edge_size=edge_size, + CGAL::parameters::facet_angle=facet_angle, + CGAL::parameters::facet_size=facet_size, + CGAL::parameters::facet_distance=facet_distance, + CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio, + CGAL::parameters::cell_size=cell_size + ); + + // Mesh generation + if (!verbose) { + // suppress output + std::cerr.setstate(std::ios_base::failbit); + } + C3t3 c3t3 = CGAL::make_mesh_3( + cgal_domain, + criteria, + lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(), + odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(), + perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(), + exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude() + ); + if (!verbose) { + std::cerr.clear(); + } + + // Output + std::ofstream medit_file(outfile); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate.hpp b/src/generate.hpp new file mode 100644 index 0000000..f49fa31 --- /dev/null +++ b/src/generate.hpp @@ -0,0 +1,33 @@ +#ifndef GENERATE_HPP +#define GENERATE_HPP + +#include "domain.hpp" + +#include +#include +#include + +namespace pygalmesh { + +void generate_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const std::vector>> & feature_edges = {}, + const double bounding_sphere_radius = 0.0, + const double boundary_precision = 1.0e-4, + const bool lloyd = false, + const bool odt = false, + const bool perturb = true, + const bool exude = true, + const double edge_size = 0.0, // std::numeric_limits::max(), + const double facet_angle = 0.0, + const double facet_size = 0.0, + const double facet_distance = 0.0, + const double cell_radius_edge_ratio = 0.0, + const double cell_size = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_HPP diff --git a/src/generate_from_off.cpp b/src/generate_from_off.cpp new file mode 100644 index 0000000..83ab7e9 --- /dev/null +++ b/src/generate_from_off.cpp @@ -0,0 +1,99 @@ +#include "generate_from_off.hpp" + +#include +#include +#include +#include +#include +#include +#include + +// IO +#include + +namespace pygalmesh { + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Polyhedral_mesh_domain_3 Mesh_domain; + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +void +generate_from_off( + const std::string & infile, + const std::string & outfile, + const bool lloyd, + const bool odt, + const bool perturb, + const bool exude, + const double edge_size, + const double facet_angle, + const double facet_size, + const double facet_distance, + const double cell_radius_edge_ratio, + const double cell_size, + const bool verbose + ) +{ + // Create input polyhedron + Polyhedron polyhedron; + std::ifstream input(infile); + input >> polyhedron; + if (!input.good()) { + std::stringstream msg; + msg << "Cannot read input file \"" << infile << "\"" << std::endl; + throw std::runtime_error(msg.str()); + } + input.close(); + + // Create domain + Mesh_domain cgal_domain(polyhedron); + + // Mesh criteria + Mesh_criteria criteria( + CGAL::parameters::edge_size=edge_size, + CGAL::parameters::facet_angle=facet_angle, + CGAL::parameters::facet_size=facet_size, + CGAL::parameters::facet_distance=facet_distance, + CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio, + CGAL::parameters::cell_size=cell_size + ); + + // Mesh generation + if (!verbose) { + // suppress output + std::cerr.setstate(std::ios_base::failbit); + } + // TODO make_mesh_3 segfaults on travis. hm. + C3t3 c3t3 = CGAL::make_mesh_3( + cgal_domain, + criteria, + lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(), + odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(), + perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(), + exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude() + ); + if (!verbose) { + std::cerr.clear(); + } + + // Output + std::ofstream medit_file(outfile); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate_from_off.hpp b/src/generate_from_off.hpp new file mode 100644 index 0000000..c2aba9c --- /dev/null +++ b/src/generate_from_off.hpp @@ -0,0 +1,28 @@ +#ifndef GENERATE_FROM_OFF_HPP +#define GENERATE_FROM_OFF_HPP + +#include +#include + +namespace pygalmesh { + +void +generate_from_off( + const std::string & infile, + const std::string & outfile, + const bool lloyd = false, + const bool odt = false, + const bool perturb = true, + const bool exude = true, + const double edge_size = 0.0, // std::numeric_limits::max(), + const double facet_angle = 0.0, + const double facet_size = 0.0, + const double facet_distance = 0.0, + const double cell_radius_edge_ratio = 0.0, + const double cell_size = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_FROM_OFF_HPP diff --git a/src/generate_surface_mesh.cpp b/src/generate_surface_mesh.cpp new file mode 100644 index 0000000..8a38b02 --- /dev/null +++ b/src/generate_surface_mesh.cpp @@ -0,0 +1,99 @@ +#define CGAL_SURFACE_MESHER_VERBOSE 1 + +#include "generate_surface_mesh.hpp" + +#include +#include +#include +#include +#include +#include + +namespace pygalmesh { + +// default triangulation for Surface_mesher +typedef CGAL::Surface_mesh_default_triangulation_3 Tr; +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2t3; +typedef Tr::Geom_traits GT; + +// Wrapper for DomainBase for translating to GT. +class CgalDomainWrapper +{ + public: + explicit CgalDomainWrapper(const std::shared_ptr & domain): + domain_(domain) + { + } + + virtual ~CgalDomainWrapper() = default; + + virtual + GT::FT + operator()(GT::Point_3 p) const + { + return domain_->eval({p.x(), p.y(), p.z()}); + } + + private: + const std::shared_ptr domain_; +}; + +typedef CGAL::Implicit_surface_3 Surface_3; + +void +generate_surface_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const double bounding_sphere_radius, + const double angle_bound, + const double radius_bound, + const double distance_bound, + const bool verbose + ) +{ + const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ? + bounding_sphere_radius*bounding_sphere_radius : + // add a little wiggle room + 1.01 * domain->get_bounding_sphere_squared_radius(); + + Tr tr; // 3D-Delaunay triangulation + C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation + + const auto d = CgalDomainWrapper(domain); + Surface_3 surface( + d, + GT::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2) + ); + + CGAL::Surface_mesh_default_criteria_3 criteria( + angle_bound, + radius_bound, + distance_bound + ); + + if (!verbose) { + // suppress output + std::cout.setstate(std::ios_base::failbit); + std::cerr.setstate(std::ios_base::failbit); + } + CGAL::make_surface_mesh( + c2t3, + surface, + criteria, + CGAL::Non_manifold_tag() + ); + if (!verbose) { + std::cout.clear(); + std::cerr.clear(); + } + + // Output + std::ofstream off_file(outfile); + CGAL::output_surface_facets_to_off(off_file, c2t3); + off_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate_surface_mesh.hpp b/src/generate_surface_mesh.hpp new file mode 100644 index 0000000..8a63043 --- /dev/null +++ b/src/generate_surface_mesh.hpp @@ -0,0 +1,23 @@ +#ifndef GENERATE_SURFACE_MESH_HPP +#define GENERATE_SURFACE_MESH_HPP + +#include "domain.hpp" + +#include +#include + +namespace pygalmesh { + +void generate_surface_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const double bounding_sphere_radius = 0.0, + const double angle_bound = 0.0, + const double radius_bound = 0.0, + const double distance_bound = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_SURFACE_MESH_HPP diff --git a/src/polygon2d.hpp b/src/polygon2d.hpp new file mode 100644 index 0000000..6a5e231 --- /dev/null +++ b/src/polygon2d.hpp @@ -0,0 +1,308 @@ +#ifndef POLYGON2D_HPP +#define POLYGON2D_HPP + +#include "domain.hpp" + +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +namespace pygalmesh { + +class Polygon2D { + public: + explicit Polygon2D(const std::vector> & _points): + points(vector_to_cgal_points(_points)) + { + } + + virtual ~Polygon2D() = default; + + std::vector + vector_to_cgal_points(const std::vector> & _points) const + { + std::vector points2(_points.size()); + for (size_t i = 0; i < _points.size(); i++) { + assert(_points[i].size() == 2); + points2[i] = K::Point_2(_points[i][0], _points[i][1]); + } + return points2; + } + + bool + is_inside(const std::array & point) + { + K::Point_2 pt(point[0], point[1]); + switch(CGAL::bounded_side_2(this->points.begin(), this->points.end(), pt, K())) { + case CGAL::ON_BOUNDED_SIDE: + return true; + case CGAL::ON_BOUNDARY: + return true; + case CGAL::ON_UNBOUNDED_SIDE: + return false; + default: + return false; + } + return false; + } + + public: + const std::vector points; +}; + + +class Extrude: public pygalmesh::DomainBase { + public: + Extrude( + const std::shared_ptr & poly, + const std::array & direction, + const double alpha = 0.0, + const double edge_size = 0.0 + ): + poly_(poly), + direction_(direction), + alpha_(alpha), + edge_size_(edge_size) + { + } + + virtual ~Extrude() = default; + + virtual + double + eval(const std::array & x) const + { + if (x[2] < 0.0 || x[2] > direction_[2]) { + return 1.0; + } + + const double beta = x[2] / direction_[2]; + + std::array x2 = { + x[0] - beta * direction_[0], + x[1] - beta * direction_[1] + }; + + if (alpha_ != 0.0) { + std::array x3; + // turn by -beta*alpha + const double sinAlpha = sin(beta*alpha_); + const double cosAlpha = cos(beta*alpha_); + x3[0] = cosAlpha * x2[0] + sinAlpha * x2[1]; + x3[1] = -sinAlpha * x2[0] + cosAlpha * x2[1]; + x2 = x3; + } + + return poly_->is_inside(x2) ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & pt: poly_->points) { + // bottom polygon + const double nrm0 = pt.x()*pt.x() + pt.y()*pt.y(); + if (nrm0 > max) { + max = nrm0; + } + + // TODO rotation + + // top polygon + const double x = pt.x() + direction_[0]; + const double y = pt.y() + direction_[1]; + const double z = direction_[2]; + const double nrm1 = x*x + y*y + z*z; + if (nrm1 > max) { + max = nrm1; + } + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features = {}; + + size_t n; + + // bottom polygon + n = poly_->points.size(); + for (size_t i=0; i < n-1; i++) { + features.push_back({ + {poly_->points[i].x(), poly_->points[i].y(), 0.0}, + {poly_->points[i+1].x(), poly_->points[i+1].y(), 0.0} + }); + } + features.push_back({ + {poly_->points[n-1].x(), poly_->points[n-1].y(), 0.0}, + {poly_->points[0].x(), poly_->points[0].y(), 0.0} + }); + + // top polygon, R*x + d + n = poly_->points.size(); + const double sinAlpha = sin(alpha_); + const double cosAlpha = cos(alpha_); + for (size_t i=0; i < n-1; i++) { + features.push_back({ + { + cosAlpha * poly_->points[i].x() - sinAlpha * poly_->points[i].y() + direction_[0], + sinAlpha * poly_->points[i].x() + cosAlpha * poly_->points[i].y() + direction_[1], + direction_[2] + }, + { + cosAlpha * poly_->points[i+1].x() - sinAlpha * poly_->points[i+1].y() + direction_[0], + sinAlpha * poly_->points[i+1].x() + cosAlpha * poly_->points[i+1].y() + direction_[1], + direction_[2] + } + }); + } + features.push_back({ + { + cosAlpha * poly_->points[n-1].x() - sinAlpha * poly_->points[n-1].y() + direction_[0], + sinAlpha * poly_->points[n-1].x() + cosAlpha * poly_->points[n-1].y() + direction_[1], + direction_[2] + }, + { + cosAlpha * poly_->points[0].x() - sinAlpha * poly_->points[0].y() + direction_[0], + sinAlpha * poly_->points[0].x() + cosAlpha * poly_->points[0].y() + direction_[1], + direction_[2] + } + }); + + // features connecting the top and bottom + if (alpha_ == 0) { + for (const auto & pt: poly_->points) { + std::vector> line = { + {pt.x(), pt.y(), 0.0}, + {pt.x() + direction_[0], pt.y() + direction_[1], direction_[2]} + }; + features.push_back(line); + } + } else { + // Alright, we need to chop the lines on which the polygon corners are + // sitting into pieces. How long? About edge_size. For the starting point + // (x0, y0, z0) height h and angle alpha, the lines are given by + // + // f(beta) = ( + // cos(alpha*beta) x0 - sin(alpha*beta) y0, + // sin(alpha*beta) x0 + cos(alpha*beta) y0, + // z0 + beta * h + // ) + // + // with beta in [0, 1]. The length from beta0 till beta1 is then + // + // l = sqrt(alpha^2 (x0^2 + y0^2) + h^2) * (beta1 - beta0). + // + const double height = direction_[2]; + for (const auto & pt: poly_->points) { + const double l = sqrt(alpha_*alpha_ * (pt.x()*pt.x() + pt.y()*pt.y()) + height*height); + assert(edge_size_ > 0.0); + const size_t n = int(l / edge_size_ - 0.5) + 1; + std::vector> line = { + {pt.x(), pt.y(), 0.0}, + }; + for (size_t i=0; i < n; i++) { + const double beta = double(i+1) / n; + const double sinAB = sin(alpha_*beta); + const double cosAB = cos(alpha_*beta); + line.push_back({ + cosAB * pt.x() - sinAB * pt.y(), + sinAB * pt.x() + cosAB * pt.y(), + beta * height + }); + } + features.push_back(line); + } + } + + return features; + }; + + private: + const std::shared_ptr poly_; + const std::array direction_; + const double alpha_; + const double edge_size_; +}; + + +class ring_extrude: public pygalmesh::DomainBase { + public: + ring_extrude( + const std::shared_ptr & poly, + const double edge_size + ): + poly_(poly), + edge_size_(edge_size) + { + assert(edge_size > 0.0); + } + + virtual ~ring_extrude() = default; + + virtual + double + eval(const std::array & x) const + { + const double r = sqrt(x[0]*x[0] + x[1]*x[1]); + const double z = x[2]; + + return poly_->is_inside({r, z}) ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & pt: poly_->points) { + const double nrm1 = pt.x()*pt.x() + pt.y()*pt.y(); + if (nrm1 > max) { + max = nrm1; + } + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features = {}; + + for (const auto & pt: poly_->points) { + const double r = pt.x(); + const double circ = 2 * 3.14159265359 * r; + const size_t n = int(circ / edge_size_ - 0.5) + 1; + std::vector> line; + for (size_t i=0; i < n; i++) { + const double alpha = (2 * 3.14159265359 * i) / n; + line.push_back({ + r * cos(alpha), + r * sin(alpha), + pt.y() + }); + } + line.push_back(line.front()); + features.push_back(line); + } + return features; + } + + private: + const std::shared_ptr poly_; + const double edge_size_; +}; + +} // namespace pygalmesh + +#endif // POLYGON2D_HPP diff --git a/src/primitives.hpp b/src/primitives.hpp new file mode 100644 index 0000000..f3e66c3 --- /dev/null +++ b/src/primitives.hpp @@ -0,0 +1,453 @@ +#ifndef PRIMITIVES_HPP +#define PRIMITIVES_HPP + +#include "domain.hpp" + +#include +#include + +namespace pygalmesh { + +class Ball: public pygalmesh::DomainBase +{ + public: + Ball( + const std::array & x0, + const double radius + ): + x0_(x0), + radius_(radius) + { + assert(x0_.size() == 3); + } + + virtual ~Ball() = default; + + virtual + double + eval(const std::array & x) const + { + const double xx0 = x[0] - x0_[0]; + const double yy0 = x[1] - x0_[1]; + const double zz0 = x[2] - x0_[2]; + return xx0*xx0 + yy0*yy0 + zz0*zz0 - radius_*radius_; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double x0_nrm = sqrt(x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2]); + return (x0_nrm + radius_) * (x0_nrm + radius_); + } + + private: + const std::array x0_; + const double radius_; +}; + + +class Cuboid: public pygalmesh::DomainBase +{ + public: + Cuboid( + const std::array & x0, + const std::array & x1 + ): + x0_(x0), + x1_(x1) + { + } + + virtual ~Cuboid() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO differentiable expression? + return std::max(std::max( + (x[0] - x0_[0]) * (x[0] - x1_[0]), + (x[1] - x0_[1]) * (x[1] - x1_[1]) + ), + (x[2] - x0_[2]) * (x[2] - x1_[2]) + ); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double x0_nrm2 = x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2]; + const double x1_nrm2 = x1_[0]*x1_[0] + x1_[1]*x1_[1] + x1_[2]*x1_[2]; + return std::max({x0_nrm2, x1_nrm2}); + } + + virtual + std::vector>> + get_features() const + { + std::vector> corners = { + {x0_[0], x0_[1], x0_[2]}, + {x1_[0], x0_[1], x0_[2]}, + {x0_[0], x1_[1], x0_[2]}, + {x0_[0], x0_[1], x1_[2]}, + {x1_[0], x1_[1], x0_[2]}, + {x1_[0], x0_[1], x1_[2]}, + {x0_[0], x1_[1], x1_[2]}, + {x1_[0], x1_[1], x1_[2]} + }; + return { + {corners[0], corners[1]}, + {corners[0], corners[2]}, + {corners[0], corners[3]}, + {corners[1], corners[4]}, + {corners[1], corners[5]}, + {corners[2], corners[4]}, + {corners[2], corners[6]}, + {corners[3], corners[5]}, + {corners[3], corners[6]}, + {corners[4], corners[7]}, + {corners[5], corners[7]}, + {corners[6], corners[7]} + }; + }; + + private: + const std::array x0_; + const std::array x1_; +}; + + +class Ellipsoid: public pygalmesh::DomainBase +{ + public: + Ellipsoid( + const std::array & x0, + const double a0, + const double a1, + const double a2 + ): + x0_(x0), + a0_2_(a0*a0), + a1_2_(a1*a1), + a2_2_(a2*a1) + { + } + + virtual ~Ellipsoid() = default; + + virtual + double + eval(const std::array & x) const + { + const double xx0 = x[0] - x0_[0]; + const double yy0 = x[1] - x0_[1]; + const double zz0 = x[2] - x0_[2]; + return xx0*xx0/a0_2_ + yy0*yy0/a1_2_ + zz0*zz0/a2_2_ - 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return std::max({a0_2_, a1_2_, a2_2_}); + } + + private: + const std::array x0_; + const double a0_2_; + const double a1_2_; + const double a2_2_; +}; + + +class Cylinder: public pygalmesh::DomainBase +{ + public: + Cylinder( + const double z0, + const double z1, + const double radius, + const double feature_edge_length + ): + z0_(z0), + z1_(z1), + radius_(radius), + feature_edge_length_(feature_edge_length) + { + assert(z1_ > z0_); + } + + virtual ~Cylinder() = default; + + virtual + double + eval(const std::array & x) const + { + return (z0_ < x[2] && x[2] < z1_) ? + x[0]*x[0] + x[1]*x[1] - radius_*radius_ : + 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double zmax = std::max({abs(z0_), abs(z1_)}); + return zmax*zmax + radius_*radius_; + } + + virtual + std::vector>> + get_features() const + { + const double pi = 3.1415926535897932384; + const size_t n = 2 * pi * radius_ / feature_edge_length_; + std::vector> circ0(n+1); + std::vector> circ1(n+1); + for (size_t i=0; i < n; i++) { + const double c = radius_ * cos((2*pi * i) / n); + const double s = radius_ * sin((2*pi * i) / n); + circ0[i] = {c, s, z0_}; + circ1[i] = {c, s, z1_}; + } + // close the circles + circ0[n] = circ0[0]; + circ1[n] = circ1[0]; + return {circ0, circ1}; + }; + + private: + const double z0_; + const double z1_; + const double radius_; + const double feature_edge_length_; +}; + + +class Cone: public pygalmesh::DomainBase +{ + public: + Cone( + const double radius, + const double height, + const double feature_edge_length + ): + radius_(radius), + height_(height), + feature_edge_length_(feature_edge_length) + { + assert(radius_ > 0.0); + assert(height_ > 0.0); + } + + virtual ~Cone() = default; + + virtual + double + eval(const std::array & x) const + { + const double rad = radius_ * (1.0 - x[2] / height_); + + return (0.0 < x[2] && x[2] < height_) ? + x[0]*x[0] + x[1]*x[1] - rad*rad : + 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double max = std::max({radius_, height_}); + return max*max; + } + + virtual + std::vector>> + get_features() const + { + const double pi = 3.1415926535897932384; + const size_t n = 2 * pi * radius_ / feature_edge_length_; + std::vector> circ0(n+1); + for (size_t i=0; i < n; i++) { + const double c = radius_ * cos((2*pi * i) / n); + const double s = radius_ * sin((2*pi * i) / n); + circ0[i] = {c, s, 0.0}; + } + circ0[n] = circ0[0]; + return {circ0}; + }; + + private: + const double radius_; + const double height_; + const double feature_edge_length_; +}; + + +class Tetrahedron: public pygalmesh::DomainBase +{ + public: + Tetrahedron( + const std::array & x0, + const std::array & x1, + const std::array & x2, + const std::array & x3 + ): + x0_(Eigen::Vector3d(x0[0], x0[1], x0[2])), + x1_(Eigen::Vector3d(x1[0], x1[1], x1[2])), + x2_(Eigen::Vector3d(x2[0], x2[1], x2[2])), + x3_(Eigen::Vector3d(x3[0], x3[1], x3[2])) + { + } + + virtual ~Tetrahedron() = default; + + bool isOnSameSide( + const Eigen::Vector3d & v0, + const Eigen::Vector3d & v1, + const Eigen::Vector3d & v2, + const Eigen::Vector3d & v3, + const Eigen::Vector3d & p + ) const + { + const auto normal = (v1 - v0).cross(v2 - v0); + const double dot_v3 = normal.dot(v3 - v0); + const double dot_p = normal.dot(p - v0); + return ( + (dot_v3 > 0 && dot_p > 0) || (dot_v3 < 0 && dot_p < 0) + ); + } + + virtual + double + eval(const std::array & x) const + { + // TODO continuous expression + Eigen::Vector3d pvec(x.data()); + const bool a = + isOnSameSide(x0_, x1_, x2_, x3_, pvec) && + isOnSameSide(x1_, x2_, x3_, x0_, pvec) && + isOnSameSide(x2_, x3_, x0_, x1_, pvec) && + isOnSameSide(x3_, x0_, x1_, x2_, pvec); + return a ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return std::max({ + x0_.dot(x0_), + x1_.dot(x1_), + x2_.dot(x2_), + x3_.dot(x3_) + }); + } + + virtual + std::vector>> + get_features() const + { + std::vector> pts = { + {x0_[0], x0_[1], x0_[2]}, + {x1_[0], x1_[1], x1_[2]}, + {x2_[0], x2_[1], x2_[2]}, + {x3_[0], x3_[1], x3_[2]} + }; + return { + {pts[0], pts[1]}, + {pts[0], pts[2]}, + {pts[0], pts[3]}, + {pts[1], pts[2]}, + {pts[1], pts[3]}, + {pts[2], pts[3]} + }; + }; + + private: + const Eigen::Vector3d x0_; + const Eigen::Vector3d x1_; + const Eigen::Vector3d x2_; + const Eigen::Vector3d x3_; +}; + + +class Torus: public pygalmesh::DomainBase +{ + public: + Torus( + const double major_radius, + const double minor_radius + ): + major_radius_(major_radius), + minor_radius_(minor_radius) + { + } + + virtual ~Torus() = default; + + virtual + double + eval(const std::array & x) const + { + const double r = sqrt(x[0]*x[0] + x[1]*x[1]); + return ( + (r - major_radius_)*(r - major_radius_) + x[2]*x[2] + - minor_radius_*minor_radius_ + ); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return (major_radius_ + minor_radius_)*(major_radius_ + minor_radius_); + } + + private: + const double major_radius_; + const double minor_radius_; +}; + + +class HalfSpace: public pygalmesh::DomainBase +{ + public: + HalfSpace( + const std::array & n, + const double alpha, + const double bounding_sphere_squared_radius + ): + n_(n), + alpha_(alpha), + bounding_sphere_squared_radius_(bounding_sphere_squared_radius) + { + } + + virtual ~HalfSpace() = default; + + virtual + double + eval(const std::array & x) const + { + return n_[0]*x[0] + n_[1]*x[1] + n_[2]*x[2] - alpha_; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return bounding_sphere_squared_radius_; + } + + private: + const std::array n_; + const double alpha_; + const double bounding_sphere_squared_radius_; +}; + +} // namespace pygalmesh + +#endif // PRIMITIVES_HPP diff --git a/src/pybind11.cpp b/src/pybind11.cpp new file mode 100644 index 0000000..0ade4be --- /dev/null +++ b/src/pybind11.cpp @@ -0,0 +1,265 @@ +#include "domain.hpp" +#include "generate.hpp" +#include "generate_from_off.hpp" +#include "generate_surface_mesh.hpp" +#include "polygon2d.hpp" +#include "primitives.hpp" + +#include +#include + +namespace py = pybind11; +using namespace pygalmesh; + + +// https://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python +class PyDomainBase: public DomainBase { +public: + using DomainBase::DomainBase; + + double + eval(const std::array & x) const override { + PYBIND11_OVERLOAD_PURE(double, DomainBase, eval, x); + } + + double + get_bounding_sphere_squared_radius() const override { + PYBIND11_OVERLOAD_PURE(double, DomainBase, get_bounding_sphere_squared_radius); + } + + // std::vector>> + // get_features() const override { + // PYBIND11_OVERLOAD( + // std::vector>>, + // DomainBase, + // get_features, + // 0.0 + // ); + // } +}; + + +PYBIND11_MODULE(_pygalmesh, m) { + // m.doc() = "documentation string"; + + // Domain base. + // shared_ptr b/c of + // + py::class_>(m, "DomainBase") + .def(py::init<>()) + .def("eval", &DomainBase::eval) + .def("get_bounding_sphere_squared_radius", &DomainBase::get_bounding_sphere_squared_radius) + .def("get_features", &DomainBase::get_features); + + // Domain transformations + py::class_>(m, "Translate") + .def(py::init< + const std::shared_ptr &, + const std::array & + >()) + .def("eval", &Translate::eval) + .def("translate_features", &Translate::translate_features) + .def("get_bounding_sphere_squared_radius", &Translate::get_bounding_sphere_squared_radius) + .def("get_features", &Translate::get_features); + + py::class_>(m, "Rotate") + .def(py::init< + const std::shared_ptr &, + const std::array &, + const double + >()) + .def("eval", &Rotate::eval) + .def("rotate", &Rotate::rotate) + .def("rotate_features", &Rotate::rotate_features) + .def("get_bounding_sphere_squared_radius", &Rotate::get_bounding_sphere_squared_radius) + .def("get_features", &Rotate::get_features); + + py::class_>(m, "Scale") + .def(py::init< + std::shared_ptr &, + const double + >()) + .def("eval", &Scale::eval) + .def("scale_features", &Scale::scale_features) + .def("get_bounding_sphere_squared_radius", &Scale::get_bounding_sphere_squared_radius) + .def("get_features", &Scale::get_features); + + py::class_>(m, "Stretch") + .def(py::init< + std::shared_ptr &, + const std::array & + >()) + .def("eval", &Stretch::eval) + .def("stretch_features", &Stretch::stretch_features) + .def("get_bounding_sphere_squared_radius", &Stretch::get_bounding_sphere_squared_radius) + .def("get_features", &Stretch::get_features); + + py::class_>(m, "Intersection") + .def(py::init< + std::vector> & + >()) + .def("eval", &Intersection::eval) + .def("get_bounding_sphere_squared_radius", &Intersection::get_bounding_sphere_squared_radius) + .def("get_features", &Intersection::get_features); + + py::class_>(m, "Union") + .def(py::init< + std::vector> & + >()) + .def("eval", &Union::eval) + .def("get_bounding_sphere_squared_radius", &Union::get_bounding_sphere_squared_radius) + .def("get_features", &Union::get_features); + + py::class_>(m, "Difference") + .def(py::init< + std::shared_ptr &, + std::shared_ptr & + >()) + .def("eval", &Difference::eval) + .def("get_bounding_sphere_squared_radius", &Difference::get_bounding_sphere_squared_radius) + .def("get_features", &Difference::get_features); + + // Primitives + py::class_>(m, "Ball") + .def(py::init< + const std::array &, + const double + >()) + .def("eval", &Ball::eval) + .def("get_bounding_sphere_squared_radius", &Ball::get_bounding_sphere_squared_radius); + + py::class_>(m, "Cuboid") + .def(py::init< + const std::array &, + const std::array & + >()) + .def("eval", &Cuboid::eval) + .def("get_bounding_sphere_squared_radius", &Cuboid::get_bounding_sphere_squared_radius) + .def("get_features", &Cuboid::get_features); + + py::class_>(m, "Ellipsoid") + .def(py::init< + const std::array &, + const double, + const double, + const double + >()) + .def("eval", &Ellipsoid::eval) + .def("get_bounding_sphere_squared_radius", &Ellipsoid::get_bounding_sphere_squared_radius) + .def("get_features", &Ellipsoid::get_features); + + py::class_>(m, "Cylinder") + .def(py::init< + const double, + const double, + const double, + const double + >()) + .def("eval", &Cylinder::eval) + .def("get_bounding_sphere_squared_radius", &Cylinder::get_bounding_sphere_squared_radius) + .def("get_features", &Cylinder::get_features); + + py::class_>(m, "Cone") + .def(py::init< + const double, + const double, + const double + >()) + .def("eval", &Cone::eval) + .def("get_bounding_sphere_squared_radius", &Cone::get_bounding_sphere_squared_radius) + .def("get_features", &Cone::get_features); + + py::class_>(m, "Tetrahedron") + .def(py::init< + const std::array &, + const std::array &, + const std::array &, + const std::array & + >()) + .def("eval", &Tetrahedron::eval) + .def("get_bounding_sphere_squared_radius", &Tetrahedron::get_bounding_sphere_squared_radius) + .def("get_features", &Tetrahedron::get_features); + + py::class_>(m, "Torus") + .def(py::init< + const double, + const double + >()) + .def("eval", &Torus::eval) + .def("get_bounding_sphere_squared_radius", &Torus::get_bounding_sphere_squared_radius) + .def("get_features", &Torus::get_features); + + py::class_>(m, "HalfSpace") + .def(py::init< + const std::array &, + const double, + const double + >()) + .def("eval", &HalfSpace::eval) + .def("get_bounding_sphere_squared_radius", &HalfSpace::get_bounding_sphere_squared_radius); + + // polygon2d + py::class_>(m, "Polygon2D") + .def(py::init< + const std::vector> & + >()) + .def("vector_to_cgal_points", &Polygon2D::vector_to_cgal_points) + .def("is_inside", &Polygon2D::is_inside); + + py::class_>(m, "Extrude") + .def(py::init< + const std::shared_ptr &, + const std::array &, + const double, + const double + >(), + py::arg("poly"), + py::arg("direction"), + py::arg("alpha") = 0.0, + py::arg("edge_size") = 0.0 + ) + .def("eval", &Extrude::eval) + .def("get_bounding_sphere_squared_radius", &Extrude::get_bounding_sphere_squared_radius) + .def("get_features", &Extrude::get_features); + + py::class_>(m, "RingExtrude") + .def(py::init< + const std::shared_ptr &, + const double + >()) + .def("eval", &ring_extrude::eval) + .def("get_bounding_sphere_squared_radius", &ring_extrude::get_bounding_sphere_squared_radius) + .def("get_features", &ring_extrude::get_features); + + // functions + m.def("generate_from_off", &generate_from_off); + m.def( + "generate_mesh", &generate_mesh, + py::arg("domain"), + py::arg("outfile"), + py::arg("feature_edges") = std::vector>>(), + py::arg("bounding_sphere_radius") = 0.0, + py::arg("boundary_precision") = 1.0e-4, + py::arg("lloyd") = false, + py::arg("odt") = false, + py::arg("perturb") = true, + py::arg("exude") = true, + py::arg("edge_size") = 0.0, + py::arg("facet_angle") = 0.0, + py::arg("facet_size") = 0.0, + py::arg("facet_distance") = 0.0, + py::arg("cell_radius_edge_ratio") = 0.0, + py::arg("cell_size") = 0.0, + py::arg("verbose") = true + ); + m.def( + "generate_surface_mesh", &generate_surface_mesh, + py::arg("domain"), + py::arg("outfile"), + py::arg("bounding_sphere_radius") = 0.0, + py::arg("angle_bound") = 0.0, + py::arg("radius_bound") = 0.0, + py::arg("distance_bound") = 0.0, + py::arg("verbose") = true + ); +} diff --git a/test/test_pygalmesh.py b/test/test_pygalmesh.py new file mode 100644 index 0000000..977840a --- /dev/null +++ b/test/test_pygalmesh.py @@ -0,0 +1,680 @@ +# -*- coding: utf-8 -*- +# +import numpy +import meshio + +# pylint: disable=import-error +import pygalmesh + + +def _row_dot(a, b): + # http://stackoverflow.com/a/26168677/353337 + return numpy.einsum("ij, ij->i", a, b) + + +def compute_volumes(vertices, tets): + cell_coords = vertices[tets] + + a = cell_coords[:, 1, :] - cell_coords[:, 0, :] + b = cell_coords[:, 2, :] - cell_coords[:, 0, :] + c = cell_coords[:, 3, :] - cell_coords[:, 0, :] + + # omega = + omega = _row_dot(a, numpy.cross(b, c)) + + # https://en.wikipedia.org/wiki/Tetrahedron#Volume + return abs(omega) / 6.0 + + +def compute_triangle_areas(vertices, triangles): + e0 = vertices[triangles[:, 1]] - vertices[triangles[:, 0]] + e1 = vertices[triangles[:, 2]] - vertices[triangles[:, 1]] + + e0_cross_e1 = numpy.cross(e0, e1) + areas = 0.5 * numpy.sqrt(_row_dot(e0_cross_e1, e0_cross_e1)) + + return areas + + +def test_ball(): + s = pygalmesh.Ball([0.0, 0.0, 0.0], 1.0) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 0]) + 1.0) < 0.02 + assert abs(max(mesh.points[:, 1]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 1]) + 1.0) < 0.02 + assert abs(max(mesh.points[:, 2]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 2]) + 1.0) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 4.0 / 3.0 * numpy.pi) < 0.15 + return + + +def test_balls_union(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Union([s0, s1]) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.1 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < 0.02 + assert abs(min(mesh.points[:, 0]) + (radius + displacement)) < 0.02 + assert abs(max(mesh.points[:, 1]) - radius) < 0.02 + assert abs(min(mesh.points[:, 1]) + radius) < 0.02 + assert abs(max(mesh.points[:, 2]) - radius) < 0.02 + assert abs(min(mesh.points[:, 2]) + radius) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 2 * ( + 4.0 / 3.0 * numpy.pi * radius ** 3 - h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2) + ) + + assert abs(vol - ref_vol) < 0.1 + + return + + +def test_balls_intersection(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Intersection([s0, s1]) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.1 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - (radius - displacement)) < 0.02 + assert abs(min(mesh.points[:, 0]) + (radius - displacement)) < 0.02 + assert abs(max(mesh.points[:, 1]) - a) < 0.02 + assert abs(min(mesh.points[:, 1]) + a) < 0.02 + assert abs(max(mesh.points[:, 2]) - a) < 0.02 + assert abs(min(mesh.points[:, 2]) + a) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 2 * (h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2)) + + assert abs(vol - ref_vol) < 0.1 + + return + + +# pylint: disable=too-many-locals +def test_balls_difference(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Difference(s0, s1) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.15 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + facet_angle=25, + facet_size=0.15, + cell_radius_edge_ratio=2.0, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + tol = 0.02 + assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < tol + assert abs(min(mesh.points[:, 0]) - 0.0) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - radius) < tol + assert abs(min(mesh.points[:, 2]) + radius) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 4.0 / 3.0 * numpy.pi * radius ** 3 - 2 * h * numpy.pi / 6.0 * ( + 3 * a ** 2 + h ** 2 + ) + + assert abs(vol - ref_vol) < 0.05 + + return + + +def test_cuboids_intersection(): + c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5]) + c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2]) + u = pygalmesh.Intersection([c0, c1]) + + # In CGAL, feature edges must not intersect, and that's a problem here: The + # intersection edges of the cuboids share eight points with the edges of + # the tall and skinny cuboid. + # eps = 1.0e-2 + # extra_features = [ + # [[1.0, 1.0 + eps, 0.5], [1.0, 2.0 - eps, 0.5]], + # [[1.0 + eps, 2.0, 0.5], [2.0 - eps, 2.0, 0.5]], + # [[2.0, 2.0 - eps, 0.5], [2.0, 1.0 + eps, 0.5]], + # [[2.0 - eps, 1.0, 0.5], [1.0 + eps, 1.0, 0.5]], + # ] + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + # filter the vertices that belong to cells + verts = mesh.points[numpy.unique(mesh.cells["tetra"])] + + tol = 1.0e-2 + assert abs(max(verts[:, 0]) - 2.0) < tol + assert abs(min(verts[:, 0]) - 1.0) < tol + assert abs(max(verts[:, 1]) - 2.0) < tol + assert abs(min(verts[:, 1]) - 1.0) < tol + assert abs(max(verts[:, 2]) - 0.5) < 0.05 + assert abs(min(verts[:, 2]) + 0.5) < 0.05 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1.0) < 0.05 + + return + + +def test_cuboids_union(): + c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5]) + c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2]) + u = pygalmesh.Union([c0, c1]) + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + # filter the vertices that belong to cells + verts = mesh.points[numpy.unique(mesh.cells["tetra"])] + + tol = 1.0e-2 + assert abs(max(verts[:, 0]) - 3.0) < tol + assert abs(min(verts[:, 0]) - 0.0) < tol + assert abs(max(verts[:, 1]) - 3.0) < tol + assert abs(min(verts[:, 1]) - 0.0) < tol + assert abs(max(verts[:, 2]) - 2.0) < tol + assert abs(min(verts[:, 2]) + 2.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 12.0) < 0.1 + return + + +def test_cuboid(): + s0 = pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]) + pygalmesh.generate_mesh(s0, "out.mesh", edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1.0) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +def test_cone(): + base_radius = 1.0 + height = 2.0 + edge_size = 0.1 + s0 = pygalmesh.Cone(base_radius, height, edge_size) + pygalmesh.generate_mesh( + s0, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 2.0e-1 + assert abs(max(mesh.points[:, 0]) - base_radius) < tol + assert abs(min(mesh.points[:, 0]) + base_radius) < tol + assert abs(max(mesh.points[:, 1]) - base_radius) < tol + assert abs(min(mesh.points[:, 1]) + base_radius) < tol + assert abs(max(mesh.points[:, 2]) - height) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = numpy.pi * base_radius * base_radius / 3.0 * height + assert abs(vol - ref_vol) < tol + return + + +def test_cylinder(): + radius = 1.0 + z0 = 0.0 + z1 = 1.0 + edge_length = 0.1 + s0 = pygalmesh.Cylinder(z0, z1, radius, edge_length) + pygalmesh.generate_mesh( + s0, "out.mesh", cell_size=0.1, edge_size=edge_length, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-1 + assert abs(max(mesh.points[:, 0]) - radius) < tol + assert abs(min(mesh.points[:, 0]) + radius) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - z1) < tol + assert abs(min(mesh.points[:, 2]) + z0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = numpy.pi * radius * radius * (z1 - z0) + assert abs(vol - ref_vol) < tol + return + + +def test_tetrahedron(): + s0 = pygalmesh.Tetrahedron( + [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0] + ) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-4 + assert abs(max(mesh.points[:, 0]) - 1.0) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 1.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1.0 / 6.0) < tol + return + + +def test_torus(): + major_radius = 1.0 + minor_radius = 0.5 + s0 = pygalmesh.Torus(major_radius, minor_radius) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-2 + radii_sum = major_radius + minor_radius + assert abs(max(mesh.points[:, 0]) - radii_sum) < tol + assert abs(min(mesh.points[:, 0]) + radii_sum) < tol + assert abs(max(mesh.points[:, 1]) - radii_sum) < tol + assert abs(min(mesh.points[:, 1]) + radii_sum) < tol + assert abs(max(mesh.points[:, 2]) - minor_radius) < tol + assert abs(min(mesh.points[:, 2]) + minor_radius) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = (numpy.pi * minor_radius * minor_radius) * (2 * numpy.pi * major_radius) + assert abs(vol - ref_vol) < 1.0e-1 + return + + +def test_custom_function(): + class Hyperboloid(pygalmesh.DomainBase): + def __init__(self, edge_size): + super(Hyperboloid, self).__init__() + self.z0 = -1.0 + self.z1 = 1.0 + self.waist_radius = 0.5 + self.edge_size = edge_size + return + + def eval(self, x): + if self.z0 < x[2] and x[2] < self.z1: + return x[0] ** 2 + x[1] ** 2 - (x[2] ** 2 + self.waist_radius) ** 2 + return 1.0 + + def get_bounding_sphere_squared_radius(self): + z_max = max(abs(self.z0), abs(self.z1)) + r_max = z_max ** 2 + self.waist_radius + return r_max * r_max + z_max * z_max + + def get_features(self): + radius0 = self.z0 ** 2 + self.waist_radius + n0 = int(2 * numpy.pi * radius0 / self.edge_size) + circ0 = [ + [ + radius0 * numpy.cos((2 * numpy.pi * k) / n0), + radius0 * numpy.sin((2 * numpy.pi * k) / n0), + self.z0, + ] + for k in range(n0) + ] + circ0.append(circ0[0]) + + radius1 = self.z1 ** 2 + self.waist_radius + n1 = int(2 * numpy.pi * radius1 / self.edge_size) + circ1 = [ + [ + radius1 * numpy.cos((2 * numpy.pi * k) / n1), + radius1 * numpy.sin((2 * numpy.pi * k) / n1), + self.z1, + ] + for k in range(n1) + ] + circ1.append(circ1[0]) + return [circ0, circ1] + + edge_size = 0.12 + d = Hyperboloid(edge_size) + + pygalmesh.generate_mesh( + d, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + # TODO check the reference values + tol = 1.0e-1 + assert abs(max(mesh.points[:, 0]) - 1.4) < tol + assert abs(min(mesh.points[:, 0]) + 1.4) < tol + assert abs(max(mesh.points[:, 1]) - 1.4) < tol + assert abs(min(mesh.points[:, 1]) + 1.4) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 1.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 2 * numpy.pi * 47.0 / 60.0) < 0.15 + return + + +def test_scaling(): + alpha = 1.3 + s = pygalmesh.Scale(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), alpha) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1 * alpha) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2 * alpha) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3 * alpha) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0 * alpha ** 3) < tol + return + + +def test_stretch(): + alpha = 2.0 + s = pygalmesh.Stretch(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [alpha, 0.0, 0.0]) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - alpha) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 12.0) < tol + + return + + +def test_rotation(): + s0 = pygalmesh.Rotate( + pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0], numpy.pi / 12.0 + ) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +def test_translation(): + s0 = pygalmesh.Translate(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0]) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 2.0) < tol + assert abs(min(mesh.points[:, 0]) - 1.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +# # segfaults on travis, works locally +# def test_off(): +# pygalmesh.generate_from_off( +# 'elephant.off', +# 'out.mesh', +# facet_angle=25.0, +# facet_size=0.15, +# facet_distance=0.008, +# cell_radius_edge_ratio=3.0, +# verbose=False +# ) +# +# vertices, cells, _, _, _ = meshio.read('out.mesh') +# +# tol = 1.0e-3 +# assert abs(max(vertices[:, 0]) - 0.357612477657) < tol +# assert abs(min(vertices[:, 0]) + 0.358747130015) < tol +# assert abs(max(vertices[:, 1]) - 0.496137874959) < tol +# assert abs(min(vertices[:, 1]) + 0.495301051456) < tol +# assert abs(max(vertices[:, 2]) - 0.298780230629) < tol +# assert abs(min(vertices[:, 2]) + 0.300472866512) < tol +# +# vol = sum(compute_volumes(vertices, cells['tetra'])) +# assert abs(vol - 0.044164693065) < tol +# return + + +def test_extrude(): + p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) + domain = pygalmesh.Extrude(p, [0.0, 0.3, 1.0]) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 0.5) < tol + assert abs(min(mesh.points[:, 0]) + 0.5) < tol + assert abs(max(mesh.points[:, 1]) - 0.8) < tol + assert abs(min(mesh.points[:, 1]) + 0.3) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 0.4) < tol + return + + +def test_extrude_rotate(): + p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) + edge_size = 0.1 + domain = pygalmesh.Extrude(p, [0.0, 0.0, 1.0], 0.5 * 3.14159265359, edge_size) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 0.583012701892) < tol + assert abs(min(mesh.points[:, 0]) + 0.5) < tol + assert abs(max(mesh.points[:, 1]) - 0.5) < tol + assert abs(min(mesh.points[:, 1]) + 0.583012701892) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 0.4) < 0.05 + return + + +def test_ring_extrude(): + p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]]) + edge_size = 0.1 + domain = pygalmesh.RingExtrude(p, edge_size) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1.5) < tol + assert abs(min(mesh.points[:, 0]) + 1.5) < tol + assert abs(max(mesh.points[:, 1]) - 1.5) < tol + assert abs(min(mesh.points[:, 1]) + 1.5) < tol + assert abs(max(mesh.points[:, 2]) - 0.5) < tol + assert abs(min(mesh.points[:, 2]) + 0.3) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 2 * numpy.pi * 0.4) < 0.05 + return + + +# def test_heart(): +# class Heart(pygalmesh.DomainBase): +# def __init__(self, edge_size): +# super(Heart, self).__init__() +# return +# +# def eval(self, x): +# return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \ +# - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3 +# +# def get_bounding_sphere_squared_radius(self): +# return 10.0 +# +# edge_size = 0.1 +# d = Heart(edge_size) +# +# pygalmesh.generate_mesh( +# d, +# 'out.mesh', +# cell_size=0.1, +# edge_size=edge_size, +# # odt=True, +# # lloyd=True, +# # verbose=True +# ) +# +# return + + +def test_sphere(): + radius = 1.0 + s = pygalmesh.Ball([0.0, 0.0, 0.0], radius) + pygalmesh.generate_surface_mesh( + s, + "out.off", + angle_bound=30, + radius_bound=0.1, + distance_bound=0.1, + verbose=False, + ) + + mesh = meshio.read("out.off") + + tol = 1.0e-2 + assert abs(max(mesh.points[:, 0]) - radius) < tol + assert abs(min(mesh.points[:, 0]) + radius) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - radius) < tol + assert abs(min(mesh.points[:, 2]) + radius) < tol + + areas = compute_triangle_areas(mesh.points, mesh.cells["triangle"]) + surface_area = sum(areas) + assert abs(surface_area - 4 * numpy.pi * radius ** 2) < 0.1 + return + + +def test_halfspace(): + c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1]) + s = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 1.0, 2.0) + u = pygalmesh.Intersection([c, s]) + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1 / 750) < 1.0e-3 + return + + +if __name__ == "__main__": + test_sphere()